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In this work, we analytically examine the validity of molecular dynamics for a soft potential 
system by considering a simple one-dimensional system with a piecewise continuous linear repulsive 
potential wall having a constant slope a. We derive an explicit analytical expression for an inevitable 
energy change due to the discrete process, which is dependent on two parameters: 1) a, which 
is a fraction of time step r immediately after the collision with the potential wall, and 2) /j = 
where po is the momentum immediately before the collision. The whole space of parameters a and 

can be divided into an infinite number of regions, where each region creates a positive or negative 
energy change AE. On the boundaries of these regions, energy does not change, ie, AE = 0. The 
envelope of \AE\ vs. fi shows a power law behavior \AE\ oc /x^, with the exponent /3 ~ 0.95. This 
implies that the round-off error in energy introduced by the discreteness is nearly proportional to 
the discrete time step r. 
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Hard, or completely impenetrable, core particles have 
an excluded volume effect. To simulate a hard-core sys- 
tem, we use event-driven methods [l]-[2] to determine the 
time at which the hard-core particles collide. At the time 
of collision, the longitudinal component of the velocities 
is exchanged, but the transverse component of the ve- 
locities remains the same after the collision. Therefore, 
their total energy is completely conserved, and there is 
absolutely no energy drift. 

In reality, atoms usually have a soft repulsive core. 
This can be modeled by exponential or power functional 
forms [3[-0. For example, noble gases such as argons 
can be described by the Lennard- Jones potential 
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where e and a characterize the energy and length scales, 
respectively. 

In molecular dynamics simulations [2], we usually cal- 
culate the positions and velocities of the particles by solv- 
ing the Newtonian equations of motion, 
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for all the particles in the system. The system governed 
by the Newtonian equations of motion has some inter- 
esting properties such as time reversibility, energy con- 
servation, and so on [5]. Molecular dynamics is a digital 
computing implementation of this continuous differential 
equation into its corresponding discrete difference equa- 
tion. Thus, keeping energy constant is a key measure in 
assessing the validity of molecular dynamics. 

An inevitable loss of accuracy is caused by represent- 
ing a derivative by its finite-difference approximation, 
termed a round-off error in digital computing. Many dif- 
ferent molecular dynamics algorithms have been 
proposed to solve Eq. (|2]). 



If the total energy of the system is different from the 
initial energy, the system may show very different ther- 
modynamic and dynamic behaviors, which no longer cor- 
respond to the original objective. Therefore, it is im- 
portant to know how the total energy changes with the 
system control parameters in the simulations. 

In this article, we derive an analytical expression for 
energy change AE in a very simple system. The system 
is composed of a simple particle colliding with the soft 
potential wall in one-dimension. The soft potential wall 
is modeled as a piecewise linear potential, as shown in 

Fig.m 
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0, ifq> 0, 
aq^ otherwise. 
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where g is a coordinate of the particle and a < is a 
constant that characterizes the slope of the potential. For 
simplicity, we take the mass of the particle as m = 1. In 

region I, where g > 0, the particle moves freely with 
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the Hamiltonian H{p^q) = which is quite trivial. In 
region II, where g < 0, it moves with the Hamiltonian 

H{p,q) = ^+aq. 

The velocity Verlet algorithm, a symplectic algorithm, 
is used to solve Eq. (|2]). For region I, it is written as 



qo -i 
Po, 
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and for region II, it is written as 
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A typical discrete trajectory for a = 0.1, r = 0.5, start- 
ing at g'o = 5.1,po = —1.0, is shown in Fig. [2l Note that 
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FIG. 1: (Color on-line) Potential model for a soft potential 
wall with slope a = 0.1. If a particle has an energy Eo = 
0.5, continuous dynamics predicts that the particle will reflect 
at q = —5.0, but in molecular dynamics, the particle may 
penetrate further into the soft wall. 
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FIG. 2: (Color on-line) A typical particle trajectory for a = 
0.1, T = 0.5 starting at = 5.1 with po = —1.0. The number 
n is the number of steps immediately before the collision with 
the potential wall when the particle enters from region I to 
II. The number n + ric is the number of steps immediately 
after the collision with the potential wall when the particle is 
leaving region II. This corresponds to a = 0.9. Note that in 
region I, the trajectory is a straight line, denoted with a star, 
but in region II, it is a parabola, denoted with a filled sphere. 



FIG. 3: (Color on-line) Divisions of a and fi space up to Uc — 
10. These divisions can continue up to ric ^ oo. The cross 
points with abscissa are given by the formula /j, = ^^-i ' ~ 
2, • • • , oo. All points on the vertical lines and diagonal lines, 
except a = show AE — 0. Upper right-angled triangles 
show > 0, denoted with the + symbol shown in the box. 
Lower right-angled triangles show AE < 0, denoted with the 
— symbol shown in the box. The numbers shown in the box 
are the number ric. 
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FIG. 4: (Color on-line) Energy change AE vs. parameter a 
for different values of /x. This corresponds to seeing a change 
along the vertical lines on the regions '4 -' and '5 +' in Fig. 
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in region I, the trajectory is a straight line, denoted with 
a star symbol, whereas in region II, the trajectory is a 
parabola, denoted with a filled sphere symbol. We now 
examine the beginning and end of the collision with the 
wall. As the particle enters a line q = from region I to 
region II, it experiences a constant repulsive force from 



the wall and then finally, leaves the region 11. Let n de- 
note the number of steps immediately before the collision 
and n-\-nc denote the number of steps immediately after 
the collision. Then, Uc would be the number of steps that 
occur during the collision. As shown in Fig. [2l when the 
particle enters the potential wall, it has a different value 
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FIG. 5: (Color on-line) Energy change vs. parameter 
fi along the horizontal lines in Fig. [3] for different values of 
parameter a. Upper panel: AE vs. /x. Lower panel: log \AE\ 
vs. log fi. The straight line shown in the lower panel has a 
slope of 0.95. 



of ^n+i, depending on the previous position Qn. We pa- 
rameterize this as a, the fraction of time during which 
the particle moves in region II from Qn to ^n+i- Thus, 
a = corresponds to g^+i = 0, and a = 1 corresponds 
to Qn = 0. The parameter a can take a value in (0, 1]. 
The case shown in Fig. [2] corresponds to a = 0.9. In 
the simulations, we actually encounter many different re- 
alizations of a; however, unfortunately, we do not know 
the distribution of a a priori. 

We derive analytical expressions for the number of col- 
lisions Tie and the relative energy change during the col- 
lision, defined by 



Pn-\-nc Pn 
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The number of collisions Uc is calculated from the in- 
equality 

{a - l)rpn + ricTpn - ^""^^^^ — ^ar^ > 0, (7) 
so that Ur can be obtained from 
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where /i = ^ and the function Ceilix) denotes the mini- 
mum integer larger than x. The quantity ar corresponds 
to the magnitude of momentum generated by the force 
field in a time step r, and thus, parameter /i is the ratio 
of this momentum and the momentum of the incoming 
particle. The relative energy change AE is calculated by 



AE = C(C - 2), 
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where ^ = jiijic — 1). Note that ric and AE are functions 
of only two dimensionless parameters, a and ji. 



Parameter a is in the range of < a < 1, and ji is in 
the range of < /i < 1 since Uc is an integer greater than 
1. From Eq. ([8j), for a fixed value of /i, the number of 
collisions Uc changes by 1 over the whole range of a from 
to 1, since nc(a = 1) — nc{oL = 0) = 1 for any value of 
ji. If ^ = /i(^c — 1) = 2, then AE = for any value of a. 
Thus, vertical lines /i = ^t^ty foi" = 2, • • • , oo divide 
the parameter space into an infinite number of subspaces, 
as shown in Fig. [3l Straight lines a = {ric — l){^ncfi — 1) 
for nc = 2, •• • , oo divide the subspaces into two regions: 
upper triangles show a positive change in AE^ and lower 
triangles show a negative change in A^. At points on the 
boundaries, AE = 0. In Fig. [3l the regions are denoted 
with different colors and a pair of symbols including a 
number Uc and a symbol, or up to Uc = 10. 

Energy changes AE as a function of a for different 
values of /i are illustrated in Fig. SI These correspond 
to seeing a change of AE along the vertical lines on the 
regions '4 -', and '5 +' in Fig. HJ For the case of /i = 0.4, 
every point except a = gives AE = 0. For other cases, 
AE = occurs when crossing the diagonal lines. 

Fig. [5] shows the results of AE as a function of /i for 
different values of a. The upper panel shows AE vs. fi 
and the lower panel shows log|A£^| vs. log/i. We find 
that the envelope of log|A£^| is linear with a slope of 
0.95. This exponent governs the energy fluctuation from 
the round-off error of a discrete time step r. 

Although we have demonstrated the validity of molec- 
ular dynamics from a very simple system for soft matter, 
the results can be generalized to a wide range of systems. 
Furthermore, if we recognize many realizations of colli- 
sions, then energy fluctuations, or energy drifts in the 
system are coming from distributions of just two vari- 
ables: 1) a, which is related to the phase of the collision 
with the potential wall when the particle enters it and 2) 
Pn-, which is different from collision to collision. 

In summary, we considered a very simple system com- 
posed of a particle entering a linear repulsive potential 
wall with a constant slope a. By applying the velocity 
Verlet algorithm, the incoming momentum Pn and the 
outgoing momentum pn+uc could be derived. Further, we 
also derived the collision number Uc and energy change 
AE. These two values depended on only two parame- 
ters, a and /i, with ranges of < a < 1 and < /i < 1, 
respectively. The parameter space created by a and /i 
could be divided into regions by the values Uc and signs 
of AE. An energy change of 0, i.e., AE = 0, was only 
seen on the boundaries of the regions. The envelope of 
log I A£^| followed a power law with an exponent of 0.95. 
Roughly speaking, |A£^| ^ /i^-^^. Finally, any two parti- 
cles with a soft repulsive potential core could be treated 
in the same way if we introduced the reduced coordinate 
system. This needs to be examined further. 
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